A  Technique  for  Developing  an  Entropy  Consistent 
System  of  Second-Order  Hydrodynamic  Equations 

Ramesh  Balakrishnan^  and  Ramesh  K.  Agarwal’^’^ 


t 


Center  for  Simulation  of  Advanced  Rockets,  University  of  Illinois  at  Urbana-Champaign,  Urbana,  IL  61801 
National  Institute  for  Aviation  Research,  Wichita  State  University,  Wichita,  KS  67260 


Abstract.  This  paper  presents  the  development  of  a  novel  set  of  second-order  hydrodynamic  equations  known  as  the 
BGK-Bumett  equations  for  computing  flows  in  the  continuum-transition  regime.  The  second-order  distribution  function 
that  forms  the  basis  of  this  formulation  is  approximated  by  the  first  three  terms  of  the  Chapman-Enskog  expansion.  Such 
and  expression,  however,  does  not  readily  satisfy  the  moment  closure  property.  Hence  an  exact  closed  form  expression 
for  the  same  is  obtained  by  enforcing  moment  closure  and  solving  a  system  of  algebraic  equations  to  determine  the 
closure  coefficients.  Through  a  series  of  conjectures  the  closure  coefficients  are  designed  to  move  the  resulting  system 
of  hydrod3mamic  equations  towards  an  entropy  consistent  set.  An  important  step  in  the  formulation  of  the  higher-order 
distribution  functions  is  the  proper  representation  of  the  material  derivatives  in  terms  of  the  spatial  derivatives.  While  the 
material  derivatives  in  the  first-order  distribution  function  are  approximated  by  the  Euler  equations,  proper  representations 
to  these  derivatives  in  the  second-order  distribution  function  are  determined  by  an  entropy  consistent  relaxation  technique. 
The  BGK-Bumett  equations,  obtained  by  taking  moments  of  the  Boltzmann  equation  in  the  second-order  distribution 
function,  are  shown  to  be  stable  to  small  wavelength  disturbances  and  entropy  consistent  for  a  wide  range  of  grid  points 
and  Mach  numbers. 


The  Need  for  Higher-Order  Hydrodynamic  Equations 

Orbital  Transfer  Vehicles  (OTVs)  belong  to  a  class  of  hypersonic  vehicles  that  are  required  to  return  to  a  low  earth  orbit 
from  a  high  earth  orbit  as  part  of  their  mission.  Consequently,  they  have  a  substantial  portion  of  their  flight  envelope  in 
the  continuum-transition  regime  which  lies  between  the  continuum  and  free  molecular  regimes.  In  this  regime  the  drag 
and  aerodynamic  heating  are  very  sensitive  to  the  degree  of  rarefaction  and  their  prediction  has  posed  quite  a  challenge 
to  designers  who  have  had  to  resort  to  empirical  correlations  based  on  sparse  experimental  data.  Since  the  ability  to 
conduct  ground  based  experiments  to  acquire  data  is  a  prohibitively  expensive  option,  it  would  be  highly  desirable  to 
have  a  computational  technique  which  can  provide  data  that  compares  favorably  with  flow  measurements  from  the  Space 
Shuttle.  Such  a  technique  may  then  be  used  to  predict  the  hypersonic  flowfield  for  future  applications. 

In  the  continuum-transition  flow  regime,  also  known  as  the  transitional  flow  regime,  the  Knudsen  number  (Kn  — 
A/Lref  where  A  denotes  the  mean  fee  path  and  Lref  denotes  the  reference  length)  is  in  the  neighborhood  of  unity.  In  this 
regime  the  Navier-Stokes  equations  yield  inaccurate  results,  as  the  approximations  made  in  the  constitutive  relations  for 
the  stress  and  heat  flux  terms,  while  acceptable  in  the  continuum  regime,  are  not  appropriate  in  the  transitional  regime. 
The  insufficient  number  of  collisions  between  the  molecules  prevents  the  gas  from  attaining  thermodynamic  equilibrium. 
This  gives  rise  to  regions  of  non-equilibrium  where  more  general  constitutive  relations  are  required  to  model  the  flow.  To 
complicate  matters  even  further,  there  may  be  regions  of  continuum  and  rarefaction  that  occur  side  by  side.  For  instance, 
in  the  flow  field  around  the  OTVs  when  they  re-enter  the  upper  atmosphere,  the  region  close  to  the  fore-body  can  be 
represented  as  a  continuum,  while  the  wake  region  exhibits  a  high  degree  of  rarefaction. 

Currently  the  only  viable  technique  for  computing  flows  in  this  regime  is  the  Direct  Simulation  Monte  Carlo  (DSMC), 
although  the  large  number  of  molecules  required  for  meaningful  results  makes  this  method  prohibitive  with  regard  to 
computational  time  and  storage  requirements.  Hence,  there  is  a  need  for  an  extended  set  of  governing  equations  which 
incorporates  more  general  expressions  for  the  constitutive  relations.  On  including  these  constitutive  relations  in  traditional 
CFD  solvers  (Navier-Stokes  solvers  in  the  continuum  domain)  it  is  expected  that  in  addition  to  capturing  the  intricacies  of 
the  flow  field  they  will  also  prove  to  be  computationally  faster  than  Monte  Carlo  simulations. 


CP585,  Rarefied  Gas  Dynamics:  22"^  International  Symposium,  edited  by  T.  J.  Bartel  and  M.  A.  Gallis 
©  2001  American  Institute  of  Physics  0-7354-0025-3/01/$18.00 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  074-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1 204,  Arlington, 
VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503 _ 

1.  AGENCY  USE  ONLY  (Leave  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

blank) _ _  9Jul  2000 _ 

4.  TITLE  AND  SUBTITLE  I  5.  FUNDING  NUMBERS 


A  Technique  for  Developing  an  Entropy  Consistent 
System  of  Second-Order  Hydrodynamic  Equations 


6.  AUTHOR(S) 

Ramesh  Balakrishnan  and  Ramesh  K.  Agarwal 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Illinois  at  Urbana- 
Champaign,  Urbana, 

Center  for  Simulation  of  Advanced 
Rockets, 

Urbana,  III. 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  Piiblic  Release;  Distribution  is  imlimited 

13.  ABSTRACT  (Maximum  200  Words) 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


15.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 


The  formulation  of  a  system  of  second-order  hydrodynamic  equations  relies  on  the  fact  that  these  equations  can  be 
obtained  by  taking  moments  of  the  Boltzmann  equation,  in  the  second-order  distribution  function,  with  the  collision  in¬ 
variant  vector.  In  an  attempt  to  derive  an  expression  for  the  second-order  stress,  Burnett  [1]  developed  a  method  by  which 
corrections  to  the  distribution  function  could  be  calculated  to  any  degree  of  approximation.  This  development  consid¬ 
ered  a  general  force  law  between  molecules  which  varied  inversely  as  the  power  of  their  distance.  In  a  subsequent 
paper,  Bimiett  [2]  derived  the  complete  expression  for  the  second-order  distribution  function  for  two  extreme  cases:  (a) 
Maxwellian  molecules,  for  which  the  force  law  varies  inversely  as  the  hfth  power  of  the  distance  and  (b)  molecules  which 
are  modeled  as  elastic  spheres. 

In  their  first  successful  attempt  at  computing  hypersonic  flows  using  a  second-order  set  of  governing  equations,  Fiscko 
and  Chapman  [3]  extended  the  numerical  methods  for  continuum  flow  into  the  continuum-transition  regime  by  incorporat¬ 
ing  the  Burnett  expressions  for  stress  and  heat  flux  into  standard  Navier-Stokes  solvers.  They  solved  the  hypersonic  shock 
structure  problem  by  relaxing  an  initial  solution  to  steady  state  and  obtained  solutions  for  a  monatomic  hard  sphere  gas 
and  argon.  Their  solutions,  obtained  on  coarse  grids  for  a  wide  range  of  Mach  numbers,  showed  that  the  Burnett  solutions 
were  in  close  agreement  with  Monte-Carlo  simulations  and  experimental  measurements  (see  Alsmeyer  [4]).  However, 
grid  refinement  studies  indicated  that  the  equations  became  unstable  as  the  mesh  size  was  made  progressively  finer.  This 
was  predicted  by  Bobylev  [5]  who  showed  that  the  Burnett  equations  are  unstable  to  small  wavelength  disturbances. 

In  an  effort  to  overcome  these  instabilities,  Zhong  [6]  formulated  the  augmented  Burnett  equations  by  adding  linear 
third-order  terms  from  the  super-Bumett  equations.  The  coefficients  (weights)  of  these  linear  third-order  terms  were 
determined  by  a  linearized  stability  analysis  of  the  augmented  Burnett  equations.  These  equations  were  successful  in 
computing  the  hypersonic  shock  structure  and  hypersonic  blunt  body  flows.  Their  initial  successes  were,  however,  short 
lived.  Attempts  at  computing  blunt  body  wakes  and  flat  plate  boundary  layers  with  the  augmented  Burnett  equations 
have  not  been  successful.  It  was  observed  that  these  equations  could  orient  the  flow  in  a  physically  unrealistic  manner  by 
allowing  shear  layers  to  sharpen  to  discontinuities  and  permitting  heat  flow  from  cold  to  hot  regions!  Since  this  behavior 
is  in  violation  of  the  second-law  of  thermodynamics  it  was  conjectured  that  this  entropy  inconsistency  may  indeed  be  the 
cause  of  the  computational  instability.  Following  this  line  of  thought  a  rigorous  thermodynamic  analysis  of  the  Burnett 
equations  ensued,  where  it  was  shown  by  Comeaux  [7]  that  the  Burnett  equations,  when  applied  to  the  hypersonic  shock 
structure  problem,  can  violate  the  second  law  of  thermodynamics  as  the  local  Knudsen  number  increases  above  a  critical 
limit. 


BGK-Burnett  Equations:  A  Novel  Second-Order  Formulation 

While  the  work  of  Comeaux  [7]  demonstrated  that  the  Burnett  equations  can  violate  the  second-law  of  thermodynamics, 
it  was  still  not  clear  what  factors  gave  rise  to  negative  irreversible  entropy.  One  of  the  causes  for  such  unphysical  effects 
could  arise  from  the  fact  that  the  higher-order  hydrodynamic  equations  do  not  form  a  closed  set.  Also,  as  mentioned  earlier, 
the  Burnett  coefficients  were  derived  for  Maxwellian  and  hard  sphere  gas  models.  For  a  real  gas,  the  coefficients  were 
assumed  to  lie  between  these  two  extremes  and  were  determined  using  an  interpolation  scheme  devised  by  Lumpkin  [8]. 
Further,  it  was  noticed,  by  the  first  author  (see  Balakrishnan  [9]),  that  although  the  Burnett  equations  did  take  into  account 
the  influence  of  forces  between  molecules  (appropriately  modeled)  it  did  not  incorporate  the  corresponding  higher-order 
virial  expansion  for  the  equation  of  state.  In  all  prior  attempts  at  computing  the  Burnett  equations  [6,  7],  the  ideal  gas 
law  had  been  assumed.  Another  cause  for  such  unphysical  effects  lies  in  the  improper  representation  of  the  material 
derivatives  in  the  second-order  stress  and  heat  flux  terms,  in  terms  of  the  spatial  derivatives.  Based  on  these  observations, 
the  following  objective  was  drawn  up  to  address  the  formulation  of  a  novel  set  of  second-order  hydrodynamic  equations 
that  is  designed  to  be  entropy  consistent. 

■  Objective.  In  order  to  understand  the  dynamics  of  the  process  by  which  the  solution  evolves  in  the  system  of  equa¬ 
tions  obtained  from  the  second-order  distribution  function,  it  was  decided  that  the  second-order  equations  be  formulated 
on  the  same  assumptions  used  to  derive  the  Navier-Stokes  equations.  While  such  a  formulation  does  not  consider  the  force 
of  attraction  between  molecules,  a  factor  that  contributes  largely  to  the  dynamics  of  molecular  collisions  at  low  pressures, 
it  must  be  emphasized  that  a  proof  of  entropy  consistency  at  this  fundamental  level  is  necessary  before  considering  all 
other  factors  that  characterize  flows  in  the  continuum-transition  regime. 


Boltzmann  Equation  with  the  BGK  model  for  the  Collision  Integral 

A  molecule  is  characterized  by  its  position  r,  velocity  v  and  internal  energy  e  which  together  constitutes  a  7-dimensional 
phase  space.  In  kinetic  theory,  a  gas  is  described  by  a  distribution  function  which  contains  information  regarding  the 


distribution  of  molecules,  their  velocities  and  (in  the  present  formulation)  their  internal  energy.  The  distribution  function 
expresses  the  probability  of  finding  molecules  in  the  physical  volume  d^r  whose  velocities  and  internal  energy  lie  in  the 
range  v  to  v  +  dv  and  e  to  e  +  de  respectively.  The  number  density  or  expected  number  of  such  molecules  per  unit  volume 
(in  physical  space)  is  defined  as 

n(r,<)  —  f  f  ^(t,r,v,e)(fvde  (1) 

Jr^Jr+ 

where  ^  =  ^{t,r,v,e)  :  K+  x  K®  x  x  R+  — >  IK+  defines  the  distribution  function.  The  flow  variables  that  are 
described  by  the  hydrodynamic  equations  are  averages  of  quantities  that  depend  on  the  velocity  and  internal  energy  of  the 
constituent  molecules.  These  averages  are  defined  as  moments  of  the  distribution  function. 

■  Definition.  The  moment  of  a  function,  <f)  —  v,  e)  :  1R+  x  E®  x  E+  M,  is  defined  as  the  inner  product 


dR3JR+ 


(2) 


over  the  molecular  space  S  —  { (v,  e)  | v  e  and  c  G  E+  } . 
The  average  of  0(f,  v,  e)  may  also  be  written  as 


0(r,  f)  =  e),  f(t,  r,  v,  e))  =  f  [  v,  e)f{t,  r,  v,  e)d^vde, 

Jr^Jr+ 


(3) 


where  the  normalized  distribution  function  is  defined  as  f(t,  r,  v,  c)  = 


n(T,t) 


If  the  gas  is  assumed  to  be  composed  of  a  single  species  of  identical  molecules,  then  in  the  absence  of  external  forces, 
the  equation  expressing  the  conservation  of  molecules  in  a  fixed  spatial  domain  D  C  E®  is  given  by 


£  l^n  if)  dn  +  I J,- in  if)  v)dQ  = 


coll 


On  multiplying  Eq.  (4)  throughout  by  the  molecular  mass  m,  we  obtain  for  an  infinitesimal  volume  5D  G  D, 


=  J[/(v,e),/(vi,ei)]. 


(4) 


(5) 


coll 


Owing  to  the  intractable  nature  of  the  collision  integral,  J  [/(v,  e),  /(vi,  ei)],  it  is  approximated  by  the  Bhatnagar-Gross- 
Krook  (BGK)  model 

</[/(v,e),/(vi,ei)]  =  -  /)  (6) 

where  v  denotes  the  collision  frequency,  the  local  Maxwellian  f^^^  =  p/Ig y/ ,  /  =  e  +  u^/2  +  u^/2, 
(3  =  1/(2RT),  Ux  denotes  the  x-component  of  the  fluid  velocity  (u),  Cx  =  Vx  —  Ux  is  the  x-component  of  the  thermal 
or  peculiar  velocity  (C)  and  Iq  —  {I,  denotes  the  average  internal  energy.  The  variables  that  are  conserved  in  a 
collision  process  are  expressed  by  the  collision  invariant  vector  ^  =  [1,  v,  J  +  (v  •  v)/2]^. 


Second-Order  Distribution  Function 


The  Chapman-Enskog  expansion  for  the  second-order  distribution  function  is  given  by  the  following  perturbative  expan¬ 
sion  about 

/  =  +  ■  •  ■  (7) 

where  the  local  Knudsen  number,  is  the  perturbation  parameter.  Substituting  Eq.  (7)  in  the  non-dimensional  form  of 
Eq.  (6)  and  equating  like  powers  of  the  Knudsen  number  gives 


f{i) 


ay(i-i) 

dt 


+  Vx 


ay(i-l) 

dx 


(8) 


In  formulating  higher-order  distribution  functions  one  starts  with  the  Maxwellian  distribution  function  and  obtains 
higher-order  terms  (iterates)  in  the  distribution  function  by  a  process  of  iterative  refinement  as  expressed  in  Eq.  (8). 


Since  the  Euler,  Navier-Stokes  and  other  higher-order  hydrodynamic  equations  must  have  the  same  held  vector,  Q  = 
[/9,  pu,  pe^ ,  it  follows  that  all  higher-order  terms  in  the  distribution  function,  1  .e.  Vt  >  1,  must  satisfy  the  following 

moment  closure  property. 

■  Moment  Closure  Property. 

(^,^i/{i))  =  0,V*>l.  (9) 

This  requirement  ensures  that  the  continuity  equation  remains  unchanged!  The  generic  expression  for  the  second-order 
distribution  function  is  given  by 


where 


0(1) 


i 

V 


^-icl  +  6^^cl  +  e,cl 
^0 


and  the  moment  closure  coefficients  satisfying  the  moment  closure  relation 
expressions 


('®-,/(o)0(i)(J,C'.)) 


(10) 


(11) 


0  are  given  by  the 


977  - 

“  4(9  -  77) 

237  -  2" 
2(9-77) 

^  177  - 

“  2(9  -  77) 

2 

^^3  =  -3  (3-7) 

02  —  (J‘^2  ~  1 

03  =  “(1  +  ^l) 

04  =  (3  —  7)  -H  0^3 

^5  =  -(7-  1)  -W3 

It  must  be  noted  that  there  is  no  unique  way  of  satisfying  the  requirement  of  moment  closure.  This  introduces  a  certain 
arbitrariness  in  the  formulation  of  the  second-order  distribution  function  and  gives  rise  to  a  family  of  BGK-Bumett  equa¬ 
tions.  Since  not  all  such  formulations  satisfy  the  twin  requirements  of  stability  and  entropy  consistency,  some  additional 
constraints  are  required  to  design  the  closure  coefficients  such  that  the  resulting  set  of  equations  are  entropy  consistent 
and  stable.  The  first  author,  in  his  doctoral  dissertation  has  described  the  formulation  and  properties  of  the  BGK-Bumett 
(1996),  BGK-Bumett  (1998)  and  BGK-Burnett  (1999)  equations.  The  closure  coefficients  given  in  this  paper  are  for  the 
BGK-Bumett  (1999)  equations.  In  the  formulation  of  the  BGK-Bumett  (1999)  equations  the  additional  constraint  that 
was  imposed  required  that  the  linearized  stability  plot  (see  Bobylev  [5]  and  Balakrishnan  [9])  of  the  BGK-Bumett  and 
N-S  equations  show  similar  variations  of  the  roots  of  the  characteristic  equation. 


Moments  of  the  Boltzmann  Equation 


A  mathematical  link  between  the  Boltzmann  equation  at  the  kinetic  level  and  the  hydrodynamic  equations  at  the  fluid  level 
is  established  by  taking  moments  of  the  Boltzmann  equation.  It  can  be  shown  that  moments  of  the  Boltzmann  equation 
with  the  Maxwellian  distribution  function  give  rise  to  the  Euler  equations.  Likewise,  moments  of  the  Boltzmann  equation 
with  the  first-order  distribution  function  give  rise  to  the  Navier-Stokes  equations.  On  taking  moments  of  the  Boltzmann 
equation  in  the  second-order  distribution  function,  Eq.  (10),  with  the  collision  invariant  vector  ^ 


dt 


(t;*^,/(°)+e/(i) 


=  0 


we  get  the  BGK-Burnett  equations 


dt  ^  dx  dx 


dG^ 


=  0 


(12) 


dx 


(13) 


where  the  flux  vectors  are  given  by 


G*  - 


pUx 

p  + pul 
PUx  “h  PUx^t 


0 


,  G®  - 


-Txx 

-Wx'Tx*  + 


,  G^- 


-T2 


(14) 


-UxT^x  +  ix 


The  stress  and  heat  flux  expressions  in  the  N-S  (superscript  v)  and  BGK-Bumett  (superscript  B)  fluxes  are  given  by 


a«  =  - 

dx  ’ 


7R  p  dT 


(7  —  1)  1/  ’ 


(15) 


'^{Kcr 


^  p  ^ 

P  St  dx  ^ 


( _ 

\dx) 


n,4§t^  +  !i,£ 

0^  m  dx  d 


0 


,  Q  ^  ^  ^  00  P 


(dd0\ 

\9xy 


+ 


^  d  dx 


0 


+ 


Q2  dp  dd 
d^  dx  dx 


(16) 


P_^(dd\  ff.  r,  \  P  ^d  dux 
^ d  dx  d^  dx  ^ d^  dx  ^ d^  ^^\dx)  ^  ^  d^  dx  dx 


III 


p  d'^Ux  p  dddux 

dd  dx  dx  ^  d^  dx^  *  d^  dx  dx 


IV 


H — 3 


£_^^_r)  P  dux  dv 
^  d^  dx  '*  d"^  dx  dx 


(17) 


The  presence  of  the  material  derivatives  (^/ =  d/dt  +  Uidjdxd  in  Eqs.  (16)  and  (17)  make  them  the  most  general 
expressions  for  the  BGK-Bumett  stress  and  heat  flux,  and  give  rise  to  a  large  number  of  representational  forms  depending 
on  the  approximations  used  to  express  these  derivatives  in  terms  of  spatial  derivatives.  The  material  derivatives  Sidj 
and  &Ux/^t  are  approximated  by  the  Euler  equations.  The  representation  of  terms  and  III  is  determined  by  an 


entropy  consistent  relaxation  technique.  The  inclusion  of  the  derivatives  of  the  collision  frequency  (v),  in  terms  II  and 
I,  in  the  formulation  of  the  BGK-Bumett  (1998)  and  BGK-Bumett  (1999)  equations  marks  an  important  step  in  the 


IV 


derivation.  Neglecting  these  terms  in  the  formulation  of  the  BGK-Bumett  (1996)  equations  yielded  expressions  for  the 
stress  and  heat  flux  that  had  the  same  derivatives  as  the  Burnett  equations.  Such  an  omission,  however,  leads  to  the 
creation  of  a  fictitious  viscosity  that  varies  as  a  function  of  pressure,  as  opposed  to  the  temperature,  thereby  causing  a 
viscous  imbalance  between  the  N-S  and  BGK-Bumett  flux  vectors  that  made  the  equations  stiff. 


The  Entropy  Consistent  Relaxation  Technique  and  Shock  Structure  Computations 

While  it  is  possible  to  treat  the  shock  as  a  sharp  discontinuity  and  use  the  Euler  equations  to  predict  the  macroscopic 
properties  of  the  shock  field,  the  structure  of  the  shock  per  se,  in  dilute  gases,  can  be  only  be  determined  by  solving  the 
Boltzmann  equation.  It  has  been  shown  by  Liepmann  et.  al  [10]  that  the  shock  structure  in  a  monatomic  perfect  gas 
is  contained  in  the  Boltzmann  equation.  Since  the  BGK-Bumett  equations  are  formulated  to  model  the  gas  in  a  state 
of  collisional  and  thermodynamic  non-equilibrium,  it  was  decided  to  solve  the  shock  stmcture  problem,  as  it  presents 
the  unique  possibilities  of  (a)  studying  the  behavior  of  the  BGK-Burnett  equations  by  isolating  the  effects  of  boundary 
conditions  and  (b)  arriving  at  an  entropy  consistent  approximation  for  the  material  derivatives  in  and  q^. 

The  Navier-Stokes  solution  to  the  shock  structure  problem  results  in  a  smooth  variation  of  the  flow  variables  about 
the  discontinuity  such  that  flux  equilibrium  is  restored  at  all  points  in  the  flowfield.  Since  the  Euler  and  Navier-Stokes 
equations  are  known  to  be  entropy  consistent,  the  resulting  shock  profile  does  not  contain  any  physical  anomalies.  How¬ 
ever,  for  the  BGK-Bumett  equations,  the  equilibrium  requirement  may  be  met  by  a  multitude  of  of  shock  shapes,  that  may 


include  “physically  untenable”  solutions.  A  more  meaningful  outcome  can  be  expected  by  insisting  that  the  governing 
equations  satisfy  the  second-law  of  thermodynamics  at  every  stage  of  the  solution  process  where  the  flow  evolves  from 
a  given  initial  profile.  Since  there  are  a  wide  variety  of  approximations  to  the  material  derivatives  in  the  BGK-Bumett 
fluxes,  one  must  identify  a  correct  approximation  that  accounts  for  the  differences  in  time  scales  between  the  first-  and 
second-order  fluxes  and  also  ensures  entropy  consistency.  In  order  to  determine  a  proper  approximation  for  the  material 
derivatives,  an  entropy  consistent  relaxation  technique  (ECRT)  has  been  developed  (see  Balakrishnan  [9]),  which  is  based 
on  the  following  premise: 

■  Premise.  By  considering  the  Navier-Stokes  solution  to  be  an  entropy  consistent  intermediate  solution  of  the  BGK- 
Burnett  equation  it  is  possible  to  select  approximations  to  the  material  derivatives  in  and  which  will  preserve  the 
positivity  of  the  irreversible  entropy  as  the  BGK-Burnett  solution  evolves. 

This  above  premise  is  based  on  the  observation  that  on  a  relatively  coarse  mesh,  where  the  local  Knudsen  number  is  quite 
low,  the  BGK-Bumett  solutions  are  indistinguishable  from  the  Navier-Stokes  solution.  Further,  by  considering  the  BGK- 
Bumett  solution  as  a  second-order  relaxation  of  the  Navier-Stokes  solution  it  is  entirely  justifiable  that  the  Navier-Stokes 
solution  would  have  developed  as  an  intermediate  solution  of  the  BGK-Bumett  equations  had  the  latter  been  started  on 
the  initial  profile  provided  by  the  Rankine-Hugoniot  (R-H)  relations.  Also,  tacit  in  this  assumption  is  the  restriction  that 
the  BGK-Bumett  equations  shall  at  no  time  during  the  solution  process  violate  the  second-law  of  thermodynamics.  The 
first  step  in  formulating  the  ECRT  is  to  develop  an  expression  for  the  irreversible  entropy  produced  by  the  BGK-Bumett 
equation.  An  expression  for  this,  based  on  the  Boltzmann  H-theorem  is  given  by  (see  Balakrishnan  [9]) 
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where  terms  of  0{p)  account  for  the  &  generated  by  the  N-S  terms  and  terms  of  0{p,^)  represent  the  contribution  of 
the  second-order  (BGK-Bumett)  terms.  The  E^’’^  coefficients  are  given  in  [9].  The  essence  of  the  ECRT  may  now  be 
summarized  in  the  following  question: 

■  Question.  When  viewed  as  an  intermediate  solution  of  the  BGK-Burnett  equations,  what  approximation  or  ap¬ 
proximations  of  the  material  derivatives  in  and  would  yield  an  entropy  consistent  Navier-Stokes  distribution  ? 

The  ECRT,  which  answers  the  above  question,  is  applied  in  the  following  manner:  For  a  given  free  stream  Mach  num¬ 
ber  the  shock  stmcture  is  computed  by  solving  the  Navier-Stokes  equations  in  the  given  control  volume  based  on  the  initial 
conditions  specified  by  the  R-H  relations.  Since  the  N-S  solution  may  be  considered  to  be  an  entropy  consistent  interme¬ 
diate  solution  of  the  BGK-Bumett  equations  it  is  imperative  that  this  solution  generates  a  positive  a  as  given  by  Eq.  ( 1 8). 
In  order  to  calculate  a  the  various  approximations  for  the  material  derivatives  identified  by  the  linearized  stability  analysis 
(see  Balakrishnan  [9])  are  substituted  in  Eq.  (18)  and  checked  for  positivity.  On  identifying  such  an  approximation  for 
the  material  derivative,  the  same  approximation  is  substituted  in  the  second-order  expressions  for  the  stress  and  heat  flux, 
Eqs.  (16-17)  and  the  BGK-Bumett  solution  proceeds  with  the  N-S  solution  as  the  initial  condition.  Based  on  the  ECRT 
it  was  identified  that  the  Navier-Stokes  approximation  for  {dux  j dx)  and  setting  {dT / dx)  —  0  give  rise  to 

entropy  consistent  expressions  for  the  second-order  stress  and  heat  flux.  Figures  l(a)-l(f)  show  the  solutions  obtained 
for  the  Mach  5  and  Mach  20  normal  shock  for  monatomic  argon  gas.  From  the  solutions  it  is  clear  that  the  differences 
between  the  N-S  and  BGK-Bumett  solutions  become  more  pronounced  as  the  free  stream  Mach  number  increases.  The 
same  is  observed  in  Tables  (1)  and  (2)  where  the  inverse  shock  thickness  (Ai/f)  based  on  the  density,  temperature  and 
stress  gradients  are  presented  for  the  Navier-Stokes  and  BGK-Bumett  (1999)  solutions. 

Blunt  Body  Flow  Computations 

In  order  to  consider  practical  applications,  a  system  of  2-D  BGK-Bumett  equations  was  derived  by  extending  the  entropy 
consistent  approximations  identified  in  the  1-D  formulation.  The  2-D  BGK-Burnett  equations  were  solved  on  a  coarse 
(57  X  81)  grid  for  the  hypersonic  flow  past  a  blunt  body  for  flow  conditions  representative  of  moderately  high  Knudsen 
numbers.  The  results  of  these  computations  for  a  free  stream  Mach  number  Moo  —  10  at  an  altitude  of  75  km  are  shown  in 
Figures  l(g)-l(l).  These  flow  conditions  correspond  to  a  free  stream  Knudsen  number  Knoo  =  0.1  for  a  cylinder  of  radius 


Upstream 

Mach  Number 

Inverse 

Density  Thickness 

Inverse 

Temperature  Thickness 

Inverse 

Stress  Thickness 

1.2 

6.306  X  10-2 

6.211  X  10-2 

4.825  X  10-2 

5.0 

0.576815 

0.471664 

0.308457 

10.0 

0.647143 

0.526993 

0.436082 

20.0 

0.667823 

0.542892 

0.439383 

Table  1:  Shock  thickness  for  argon  computed  using  the  Navier-Stokes  equations.  Fine  grid  solutions,  Ax/Ai  —  0.1. 


Upstream 

Inverse 

Inverse 

Inverse 

Mach  Number 

Density  Thickness 

Temperature  Thickness 

Stress  Thickness 

1.2 

6.20585  X  10-2 

6.18614  X  10-2 

9.21306  X  10-3 

5.0 

0.503101 

0.375218 

0.184188 

10.0 

0.570015 

0.405715 

0.307753 

20.0 

0.586499 

0.413692 

0.308121 

Table  2:  Shock  thickness  for  argon  computed  using  the  BGK-Bumett(  1999)  equations.  Fine  grid  solutions,  Aa;/Ai  =  0.1. 


Rb  =  0.02m.  It  is  seen  that  there  are  differences  in  the  solutions  obtained  from  the  Navier-Stokes  and  BGK-Bumett 
(1999)  in  the  region  upstream  of  the  bow  shock. 


Conclusions 

A  technique  to  develop  a  1-D  entropy  consistent  set  of  second-order  hydrodynamic  equations  has  been  presented.  These 
equations,  termed  the  BGK-Bumett  equations  have  been  shown  to  yield  thicker  shocks  when  compared  to  the  N-S  solu¬ 
tions  and  are  entropy  consistent  for  the  range  of  Mach  numbers  and  grid  densities  considered  in  this  study.  The  identi¬ 
fication  of  the  fictitious  viscosity  arising  as  a  result  of  neglecting  the  derivatives  of  v  represents  an  important  step  in  the 
formulation  of  an  entropy  consistent  set  of  equations.  This  formulation  which  is  based  on  the  assumption  that  there  is  no 
intermolecular  force  between  molecules  needs  to  be  extended  to  include  molecular  forces  that  vary  as  a  function  of  the 
intermolecular  distances.  It  is  also  shown  that  a  direct  extension  of  this  formulation  to  2-D  gives  rise  to  thicker  shocks  for 
flows  past  blunt  bodies.  However,  it  remains  to  be  seen  if  an  extension  of  this  methodology  to  higher  dimensions  results 
in  entropy  consistent  formulations. 
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Figure  1 :  Figures  (a),  (b),  and  (c)  show  the  variations  of  the  stress,  heat  flux  and  irreversible  entropy  due  to  the  Navier- 
Stokes,  and  BGK-Bumett  fluxes  and  their  total  contribution  across  a  Mach  5  normal  shock.  Figure  (d)  shows  the  tem¬ 
perature  and  density  variations  across  a  Mach  5  normal  shock  while  figures  (e)  and  (f)  show  the  same  variations  across  a 
Mach  20  normal  shock.  Figures  (g)  and  (h)  show  the  temperature  contours  for  the  Navier-Stokes  and  BGK-Bumett  equa¬ 
tions  respectively  while  figures  (i)  and  (j)  show  the  velocity  contours  for  the  Navier-Stokes  and  BGK-Bumett  equations. 
Figures  (k)  and  (1)  show  the  variations  of  the  flow  properties  along  the  stagnation  streamline. 


